! Bessel I_1(x) function in double precision
!
      function dbesi1(x)
      implicit real*8 (a - h, o - z)
      dimension a(0 : 59), b(0 : 69), c(0 : 44)
      data (a(i), i = 0, 11) / 
     &    1.2787464404046789181d-10, 3.5705860060088241077d-9, 
     &    9.9611537619347335040d-8, 2.2395070088633043177d-6, 
     &    4.0312466928887462346d-5, 5.6437387840203722356d-4, 
     &    5.9259259312934746096d-3, 4.4444444443499008870d-2, 
     &    2.2222222222232042719d-1, 6.6666666666666139867d-1, 
     &    1.0000000000000001106d0, 4.9999999999999999962d-1 / 
      data (a(i), i = 12, 23) / 
     &    1.7281952384448634449d-10, 3.0647204559976390130d-9, 
     &    1.0237662138842827028d-7, 2.2299494417341498163d-6, 
     &    4.0335364374929326943d-5, 5.6433440269141349899d-4, 
     &    5.9259754885893798654d-3, 4.4444399410880397870d-2, 
     &    2.2222225112835026730d-1, 6.6666665422146063244d-1, 
     &    1.0000000032274936821d0, 4.9999999961866867205d-1 / 
      data (a(i), i = 24, 35) / 
     &    2.3216048939948030996d-10, 1.7443372702334489579d-9, 
     &    1.1596478963485415499d-7, 2.1446755518623035147d-6, 
     &    4.0697440347437076195d-5, 5.6324394900433192204d-4, 
     &    5.9283484996093060678d-3, 4.4440673899150997921d-2, 
     &    2.2222638016852657860d-1, 6.6666358151576732094d-1, 
     &    1.0000013834029985337d0, 4.9999971643129650249d-1 / 
      data (a(i), i = 36, 47) / 
     &    3.1013758938255172562d-10, -8.4813676145611694984d-10, 
     &    1.5544980187411802596d-7, 1.7811109378708045726d-6, 
     &    4.2945322199060856985d-5, 5.5344850176852353639d-4, 
     &    5.9590327716950614802d-3, 4.4371611097707060659d-2, 
     &    2.2233578241986401111d-1, 6.6654747300463315310d-1, 
     &    1.0000756505206705927d0, 4.9997803664415994554d-1 / 
      data (a(i), i = 48, 59) / 
     &    4.1214758313965020365d-10, -5.3613317735347429440d-9, 
     &    2.4661360807517345161d-7, 6.7144593918926723203d-7, 
     &    5.1988027944945587571d-5, 5.0165568586065803067d-4, 
     &    6.1717530047005289953d-3, 4.3745229577317251404d-2, 
     &    2.2363147971477747996d-1, 6.6475469131117660240d-1, 
     &    1.0015686689447547657d0, 4.9941120439785391891d-1 / 
      data (b(i), i = 0, 13) / 
     &    6.6324787943143095845d-8, 4.5125928898466638619d-7, 
     &    6.7937793134877246623d-6, 7.4580507871505926302d-5, 
     &    7.6866382927334005919d-4, 7.1185174803491859307d-3, 
     &    5.8721838073486424416d-2, 4.2473949281714196041d-1, 
     &    2.6396965606282079123d0, 13.710008536637016903d0, 
     &    57.158647688180932003d0, 179.46182892089389037d0, 
     &    377.57997362398478619d0, 399.87313678256009819d0 / 
      data (b(i), i = 14, 27) / 
     &    1.7652713206027939711d-7, 1.1988179244834708057d-6, 
     &    1.8037851545747139231d-5, 1.9775785516370314656d-4, 
     &    2.0354870702829387283d-3, 1.8822164191032253600d-2, 
     &    1.5500485219010424263d-1, 1.1190100010560573210d0, 
     &    6.9391565185406617552d0, 35.948170579648649345d0, 
     &    149.41909525103032616d0, 467.42979492780642582d0, 
     &    979.04227423171290408d0, 1030.9147225169564443d0 / 
      data (b(i), i = 28, 41) / 
     &    4.7022299276154507603d-7, 3.1878571710170115972d-6, 
     &    4.7940153875711448496d-5, 5.2496623508411440227d-4, 
     &    5.3968661134780824779d-3, 4.9837081920693776234d-2, 
     &    4.0979593830387765545d-1, 2.9533186922862948404d0, 
     &    18.278176130722516369d0, 94.476497150189121070d0, 
     &    391.66075612645333624d0, 1221.4182034643210345d0, 
     &    2548.6177980961291004d0, 2670.9883037012546541d0 / 
      data (b(i), i = 42, 55) / 
     &    1.2535083724002034147d-6, 8.4845871420655708250d-6, 
     &    1.2753227372734042108d-4, 1.3950105363562648921d-3, 
     &    1.4325473993765291906d-2, 1.3212452778932829125d-1, 
     &    1.0849287786885151432d0, 7.8068089156260172673d0, 
     &    48.232254570679165833d0, 248.80659424902394371d0, 
     &    1029.0736929484210803d0, 3200.5629438795801652d0, 
     &    6656.7749162019607914d0, 6948.8586598121632302d0 / 
      data (b(i), i = 56, 69) / 
     &    3.3439394490599745013d-6, 2.2600596902211837757d-5, 
     &    3.3955927589987356838d-4, 3.7105306061050972474d-3, 
     &    3.8065263634919156421d-2, 3.5068223415665236079d-1, 
     &    2.8760027832105027316d0, 20.665999500843274339d0, 
     &    127.47939148516390205d0, 656.43636874254000885d0, 
     &    2709.5242837932479920d0, 8407.1174233600734871d0, 
     &    17437.146284159740233d0, 18141.348781638831600d0 / 
      data (c(i), i = 0, 8) / 
     &    -2.8849790431465382128d-15, -3.5125350943844774657d-14, 
     &    -7.4850867013707419750d-13, -1.8383904048277485153d-11, 
     &    -5.7303556446977223342d-10, -2.4449502737311496525d-8, 
     &    -1.6765373351766929724d-6, -3.2189516835265773471d-4, 
     &    5.1503226936425277377d-2 / 
      data (c(i), i = 9, 17) / 
     &    -5.8674306822281631119d-15, -9.4884898451194085565d-15, 
     &    -8.5033865136600364340d-13, -1.8142997866945285736d-11, 
     &    -5.7340238386338193949d-10, -2.4449138101742183665d-8, 
     &    -1.6765375646678855842d-6, -3.2189516826945356325d-4, 
     &    5.1503226936412017608d-2 / 
      data (c(i), i = 18, 26) / 
     &    -1.4723362506764340882d-14, 1.3945147385179042899d-13, 
     &    -1.9618041857586930923d-12, -1.3343606394065121821d-11, 
     &    -5.8649674606973244159d-10, -2.4426060539669553778d-8, 
     &    -1.6765631828366988006d-6, -3.2189515191449587253d-4, 
     &    5.1503226931820146445d-2 / 
      data (c(i), i = 27, 35) / 
     &    -5.8203519372580372987d-14, 1.2266326995309845825d-12, 
     &    -1.3921625844526453237d-11, 6.2228025878281625469d-11, 
     &    -8.8636681342142794023d-10, -2.3661241616744818608d-8, 
     &    -1.6777870960740520557d-6, -3.2189402882677074318d-4, 
     &    5.1503226479551959376d-2 / 
      data (c(i), i = 36, 44) / 
     &    -4.5801527369223291722d-14, 6.7998819697143727209d-13, 
     &    -4.1624857909290468421d-12, -3.2849009406112440998d-11, 
     &    -3.2478275690431118270d-10, -2.5739209934053714983d-8, 
     &    -1.6730566573215739195d-6, -3.2190010909008684076d-4, 
     &    5.1503229866932077150d-2 / 
      w = abs(x)
      if (w .lt. 8.5d0) then
          t = w * w * 0.0625d0
          k = 12 * (int(t))
          y = (((((((((((a(k) * t + a(k + 1)) * t + 
     &        a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + 
     &        a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + 
     &        a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + 
     &        a(k + 11)) * w
      else if (w .lt. 12.5d0) then
          k = int(w)
          t = w - k
          k = 14 * (k - 8)
          y = ((((((((((((b(k) * t + b(k + 1)) * t + 
     &        b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + 
     &        b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + 
     &        b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + 
     &        b(k + 11)) * t + b(k + 12)) * t + b(k + 13)
      else
          t = 60 / w
          k = 9 * (int(t))
          y = ((((((((c(k) * t + c(k + 1)) * t + 
     &        c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + 
     &        c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + 
     &        c(k + 8)) * sqrt(t) * exp(w)
      end if
      if (x .lt. 0) y = -y
      dbesi1 = y
      end
!
